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APPARATUS FOR ANALYSING CARDIAC EVENTS 
Technical field 

The present invention relates to an apparatus for analysing cardiac events 
5 detected in electrograms, EGMs, and to a heart stimulator provided with such an 
apparatus. 

In the following the expression cardiac event denotes the depolarization 
phase in the cardiac cycle which for atrial signals is commonly known as P wave 
and for ventricular signals as R wave or QRS complex. 

10 

Background 

In the field of devices for cardiac rhythm management (CRM), accurate 
rhythm classification is an increasingly important aspect. Pacemakers are primarily 
used to assist in bradycardia or when the electrical propagation path is blocked, 

15 whereas the primary use of implantable cardioverter defibrillators (ICD) is to ter- 
minate ventricular arrhythmia, a life-threatening condition if not acutely treated. In 
both types of devices, accurate event classification of the electrogram signal is 
needed for identifying, e.g., atrial and ventricular fibrillation in order to give appro- 
priate therapy for the detected arrhythmia. For pacemakers, this may imply 

20 changing the pacing mode in order to stabilize the ventricular rhythm during an 
episode of atrial fibrillation. An ICD responds to ventricular fibrillation by giving a 
defibrillator shock which hopefully terminates the fibrillation. 

Ever increasing demands are put on both kinds of devices to better hand- 
le their primary task as well as to manage other tasks than those originally inten- 

25 ded for. One such task may, for an implantable medical device, be to identify atrial 
flutter in order to terminate it by atrial pacing or to defibrillate atrial fibrillation. 
Although it is not a life-threatening arrhythmia, atrial fibrillation is an inconvenience 
to the patient and increases the risk for other diseases such as stroke. Atrial 
pacing may also be one way of terminating supraventricular tachycardias. An ICD 

30 specific task is to identify atrial fibrillation in order to not mistake it for ventricular 
fibrillation and the risk of giving an unnecessary, and possibly harmful, 
defibrillation shock. Another, more general, utilization is to efficiently store rhythm 
data for later analysis and evaluation, already done in modem ICD's. By collecting 
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data, better knowledge of the evolution of cardiac diseases and the functionality of 
the device can be obtained. 

Clustering represents an important task within the classification problem 
where each individual event is assigned to a cluster of events with similar features. 
5 Labelling of the clusters, i.e., associating the cluster with a specific cardiac rhythm, 
completes the classification such that the device can provide proper therapy when 
needed. However, certain constraints distinguish clustering in CRM devices from 
clustering in general. In order to give immediate therapy, it requires clustering to 
be done in real-time, thus excluding many iterative clustering algorithms such as k- 

10 means clustering and competitive learning. Various methods have recently been 
presented concerning clustering of signals from the surface electrocardiogram 
(ECG), based on, e.g., self-organizing maps or fuzzy hybrid neural networks, see 
M. Lagerholm, C. Peterson, G. Braccini, L. Edenbrabdt, and L. Sornmo, Clustering 
ECG complexes using Hermite functions and self-organizing maps", IEEE Trans. 

15 Biomed. Eng., vol. 47, pp. 838 -848, July 2000, and S. Osowski and T Linh, "ECG 
beat recognition using fuzzy hybrid neural network", IEEE Trans. Biomed. Eng., 
vol. 48, pp. 1265 -1271, November 2001 . However, most clustering algorithms 
used for ECG analysis are computationally rather complex and therefore unsuitab- 
le for implantable CRM devices. Furthermore, not much a priori morphologic infor- 

20 mation is associated with the various rhythms in the electrogram (EGM); this is in 
contrast to the more well-defined ECG. 

Previously presented work in the area of intracardiac event classification 
mainly focus on discrimination of a specific condition in order to discern, e.g., atrial 
fibrillation from other atrial tachyarrythmias, see A. Schoenwald, A. Sahakian, and 

25 S. Swiryn, "Discrimination of atrial fibrillation from regular atrial rhythms by spatial 
precision of local activation direction", IEEE Trans. Biomed. Eng., vol. 44, pp. 958 
- 963, October 1997. Other applications involve discrimination of ventricular from 
supraventricular tachycardia, see L. Koyrakh, J. Gillberg, and N. Wood, "Wavelet 
based algorithms for EGM morphology discrimination for implantable ICDs", in 

30 Proc. Of Comp. In Card. (Piscataway, NJ, USA), pp. 343 - 346, IEEE, IEEE Press, 
1999, and G.Gronefeld, B. Schulte, S. Hohnloser, H. - J. Trappe, T. Korte, C. 
Stellbrink, W. Jung, M. Meesmann, D. Bocker, D. Grosse - Meininghaus, J. Vogt, 
and J. Neuzner, u Morphology discrimination: A beat-to-beat algorithm for the 
discrimination of ventricular from supraventricular tachycardia by implantable 
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cardioverter befibrillators", J. Pacing Clin. Electrophysiol., vol. 24, pp. 1519-1524, 
October 2001. More general classification algorithms, which in turn involve training 
on individual patients, have been based on analogue neural networks or wavelet 
analysis for morphologic discrimination of arrhythmias. 
5 Future pacemakers and ICD's will include more advanced means for 

arrhythmia detection and therapy, and the purpose of the present invention is to 
propose a technique for separation of cardiac rhythms in a reliable way on the 
basis of electrogram event clustering 

1 0 Disclosure of the invention 

The above purpose is obtained by an apparatus according to claim 1 . 
Certain arrhythmias are diagnosed immediately to give proper therapy, 
while, for others, it may be sufficient to record the rhythm for data collection 
purposes. Thus, the classification problem, viz. to label the rhythms based on 
15 clusters using clinical terms, may not always be necessary to implement. 

According to advantageous embodiments of the apparatus according to 
the invention both morphologic and temporal data are considered for clustering. 
Morphologic features are efficiently extracted by use of the dyadic wavelet 
transform after which the events are grouped by a leader-follower clustering 
20 embodiment. The event detection problem, based on the same transform, is 
previously treated in Swedish patent application no. 0103562-5. 

According to another advantageous embodiment of the apparatus 
according to the invention an integrating means is provided to integrate said 
distance over a predetermined period of time. By integrating the distance over a 
25 period of time it is possible to distinguish irregular rhythms, like atrial fibrillation, 
from regular rhythms. The integral total distance in case of atrial fibrillation will be 
high whereas regular rhythms will result in a lower total distance. 

The invention also relates to a heart stimulator provided with the above 
mentioned apparatus for controlling the therapeutic stimulation depending on 
30 arrhythmia detection. 

Brief description of the drawings 

To explain the invention in greater detail embodiments of the invention, 
chosen as examples, will now be described with reference to the enclosed 
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drawings, on which figure 1 shows impulse responses of a filter bank used for 
cardiac event detection, figure 2 is a flow chart of clustering algorithm perform with 
the apparatus according to the invention, figure 3 illustrates the computational 
complexity for different clustering algorithms, figure 4 presents in (a) clustering 
5 performance in terms of probability of a correct clustering of an event, Pec and 
probability of a dominant event in a cluster P oe and in (b) the spread in the number 
of clusters, figure 5 illustrates in (a) clustering performance in terms of Pde and 
Pec as a function of a distance threshold q and in (b) the numbers of clusters as a 
function of q. figure 6 exemplifies the clustering performance for one set of para- 
10 meters, figure 7 exemplifies clustering result on a concatenated EGM for three 
cases, figure 8 is a flow chart illustrating the function in broad outline of an 
embodiment of the apparatus according to the invention, and figure 9 is a block 
diagram of an exemplifying embodiment of a heart stimulator provided with an 
apparatus according to the invention. 

15 

Description of preferred embodiments 
P wave detection and feature extraction 

A signal model assuming that the event waveform is composed of a linear 
combination of representative signals is considered. The feature extraction 
20 problem is then to estimate the individual components of the representative 

signals since each morphology will have its own linear combination. By using the 
dyadic wavelet transform, different widths of the two fundamental monophasic and 
biphasic waveforms are included in the model at a tow cost. 




25 Feature extraction 

It is assumed that the QRS waveform is composed of a linear combination 
of representative signals, 

H*[hi h p ] (2) 

where each function, h j? j= 1 r t P, is of size M x 1 . The only restriction on H is 

30 that it must have full rank. Different morphologies, s(n), are modelled by the P x 1 
coefficient vector 6(n), with the linear model 

s(n)=H9(n) (3) 
where n is a temporal variable describing when the event occurs. The observed 
signal, x(n), is assumed to be modelled by, 



x(n) =H6(n) + w(/j) (4) 
where the M x 1 noise vector w(n) is assumed to be zero-mean, white, and 
Gaussian with variance a * Consequently, x(n) is defined as 

x(n) 

x(n)= (5) 
x(n + M-1) 

implying an event duration of M samples, beginning at n. The probability density 
function x(n) for a specific realization of 9(n), p(x(n); 9(n)). is thus given by 



p(x(n)0(n)) = 



exp 



— ^(x(n)-H6(n)) T (x(n)-He(n)) 
2a 

to 



(6) 



In this model, a complete description of a QRS complex is provided by the 
10 deterministic unknown parameter vector 0(n). The absence of a QRS complex 
corresponds to the case when 8(n) is equal to 0, where 0 is the P x 1 zero vector 
In general, no a priori knowledge is available on 9(n), and therefore an estimate is 
required before detection can take place. Furthermore, only one event is assumed 
to take place within the observation interval 0 < n < N. 

15 

Filterbank representation 

The descriptive functions in H have been selected such that the following 
three aspects have been taken into particular account: 

1 . the main morphologies of the QRS complex are mono- and/or biphasic 
20 2. the broad range of QRS complex durations, and 

3. a low complexity implementation. 

The wavelet transform is particularly suitable since it is a local transform, 
i.e., it provides information about the local behaviour of a signal. One wavelet 
decomposition method which may be efficiently implemented is the dyadic wavelet 
25 transform. By careful selection of the filters, a suitable filter bank including mono- 
and biphasic impulse responses can be obtained. A symmetric lowpass filter, f(n) ( 
is used repeatedly in order to achieve proper frequency bands. This filter is 
combined with one of two filters, gt>(n) or g m (n), which together define the 
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waveform morphology (the subindices b and m denote bi- and monophasic, 
respectively). For the biphasic case, the recursion is expressed as, 
hi,b(n)= g b (n) 
h 2 j>(n)=f(n)*g b (2n) 
5 h 3 ^(n)=f(n)M(2n)*g b (4n) 

h qmax , b (n) = f(n) * ... * fa*™ ~ 3 n)* Qb (2*«- " 1 n) (7) 

in which the subindex represents the maximum (coarsest) scale. 
10 It is now possible to present an expression for H which is composed of 

one biphasic and one monophasic part, 

H = [H b H m ] (8) 

15 where the biphasic H b is defined as 



H b -[h qfnjn ,b-h qmax(b ] = 



h qm , ;b (0) h W/b (0) 
h qmjn . b (M-l)...h qmaxib (M-l)J 



(9) 



where the subindex q m!n represents the minimum scale and q m tn < <W. The 
20 monophasic matrix H m is computed in a corresponding way. The reversed order of 
the columns in H, denoted with H in (8), is introduced in order to be consistent with 
the model assumed in (4). 

In order to mimic the desired mono- and biphasic waveforms with a 
low complexity filter bank structure, short filters with small integer coefficients were 
25 used. In (7), the impulse response f(n) was chosen as a third order binomial 
function, 

F(z) = (1 + z" 1 f — 1 + 3z~ 1 + 3z~ 2 + z~ 3 (10) 



where F(z) is the Z-transform of f(n). For the biphasic filter bank, the 
30 filter g b (n) was selected as the first order difference, 



gt>(n) = [-1 1] ( 11 > 
The fitter g^n) was chosen such that a compromise between the 
requirement of a DC gain equal to zero and an approximately monophasic impulse 
response was achieved. A reduction of complexity results when g b {n) is reused, 
5 gm(n) =g b (n) * gz> (n) = [1-21] (12) 

For this particular choice of g m (n) and by using Mallat's algorithm, it is 
possible to calculate both the biphasic and the monophasic filter output from each 
scale by using f(n) once and g b (n) twice. The filter bank impulse responses are 
shown in Fig. 1 . The filter bank includes two orthogonal signal sets where the 

10 width of the signal varies within each set. In (a), h j>b (n) is shown for j = 2 ,4 

from top to bottom, and in (b), the corresponding hj, m (n) are shown. 
ML parameter estimation 

The unknown coefficient vector 9(n) can be estimated by using the 
maximum likelihood criterion according to, 

1 5 G(n) = argmax p(x(n)e(n)) (13) 

9(n) 

However, 9 (n) is only of interest for those n for which the probability of an 
event, or, equivalents, for which the likelihood ratio test function T(x(n)) is 
maximized, 

n = argmax T(x(n)) (14) 
n 

20 For this case, T (x(n)) can be shown to be [14], 

T(x(n)) = x(n) T H^ T H)" 1 H T x(n) (15) 
for the case when a * is assumed to be constant. 
The optimal estimate 6 for the detected event at n is thus expressed as 
e = argm^jp(x(n)e(n)) (16) 

25 The MLE of 9(n)is found by maximizing p(x(n);9(n),or equivalently 

minimizing the MSE, 

(x(n)-He(n^ 17) 
Derivation with respect to x(n) yields the optimum 9 , 

e(n)^ T H)" 1 H T x(n)4 T Hf 1 H T x(n) (18) 
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By using the above formulation, it is also possible to derive a generalized 
likelihood ratiodetector based on (15). 

Rate 

5 The central parameter for classification of cardiac arrhythmias is the heart 

rate. Most arrhythmias are defined in terms of heart rate, although sometimes with 
rather fuzzy limits. Consequently, rate should be considered in order to Improve 
performance. The RR interval, At, is defined as the duration between two 
consecutive events, 

10 At k =(n k ^n k _ 1 )T a (19) 

where n k and n k _i denote the occurrence times of the events, and T s denotes the 
sampling period. 

Leader-follower clustering 

15 The choice of leader-follower clustering is based on a number of features 

which makes it suitable for the present invention, viz. 
o on-line processing (non-iterative), and 

o self-learning, i.e., no a priori knowledge of the of clusters is needed. 
The starting point is the assumption that an event is present for which it 
20 should be decided whether it belongs to an already existing cluster or if a new 
cluster should be initiated. Since, no knowledge is a priori available on which 
rhythms or morphologies to be expected, the chosen algorithm must be self- 
learning. The leader-follower clustering algorithm is constituted by four quantities: 

1. The event parameter vector, <D(n k ), containing the features of the k:th 
25 event that occurs at time n k . 

2. The cluster center, \if , and the covariance matrix Si that together define 
the i:th cluster. Since both n, and E,are unknown, they are replaced by their 
estimates^ and irrespectively. 

3. The metric, dj 2 , determine the distance between 0(n k ) and \x f 
30 according to some suitable function. 

4. A rule for adaptation of the cluster parameters for the winning cluster. 
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Initialization of new clusters 

During run-time, a finite number of clusters exists which represent the 
rhythms having appeared until present time. Thus, it is occasionally necessary to 
initialize a new cluster when the existing ones do not sufficiently well fit the present 

5 event. When the distance function dj 2 exceeds a certain threshold, tj, it is more 
likely that the event belongs to a new cluster than to any of the existing clusters. 
The selection of n is a tradeoff between cluster size and cluster resolution in the 
sense that choosing a small ti will result in many clusters with few clustering 
errors. On the other hand, a large n results in few clusters but in more errors. 

10 The minimum distance between O(n) and \n with respect to both i and 

n is 



d min=^n d ?(0 i = l 1 '^k-f n k +| (20) 



15 over the search interval K. The corresponding minimum distance 

indices are found as 

[imin»nmin] = argmindj 2 (l) (21) 

The decision rule based on the comparison of d mi * and r\\s expressed 

as r — 
20 >n : Initialize new cluster 

<r|: Assign to winning cluster 

<pT|: Update wining cluster (22) 



d . 2 ^< 

mm *\ 



where 0 < p <1 . When the upper relation in (22) holds, a new cluster is 
25 initialized by first increasing the number of clusters by one, 1 = 1 + 1, given that 
I < Lax where Imaxis the maximum number of clusters, and then initializing the new 
cluster as, 

Mi=«>(n k ) (23) 
On the other hand, if I = Uax the algorithm needs to discard one of the 
30 existing clusters. This can be done by elimination of, e.g., the oldest or the 
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"smallest" cluster. By the term "smallest" cluster is meant that cluster which most 
rarely is fitted with a detected cardiac event. 

For the middle and lower conditions in (22), the existing minimum distance 
cluster /min is selected as the winning cluster. However, only the lower condition 
5 results in a cluster parameter update. The reason for including such a distinction is 
that only closely similar events should be used for cluster updates in order to 
reduce contamination. 

Mahalanobls distance function 

10 The distance between the feature vector 6(n) and each cluster £ , is de- 

fined as the Mahalanobis distance, which is a normalized Euclidean distance in 
the sense that it projects the parameter vector elements onto univariate dimen- 
sions by including the inverse covariance matrix If 1 . Thus, a feature with a larger 
variance in <t>(n) will be assigned a larger share of the hyperspace before norma- 
ls lization compared to that with a lower variance. A consequence of normalization is 
that the Mahalanobis distance works welt on correlated data since If 1 then acts as 
a decorrelator. 

When searching for the minimum distance, a grid search over n is per- 
formed. This grid search is necessary since it not only minimizes the distance but 
20 also results in a more accurate fiducial point estimate than what would be the case 
when only considering T(x(n k )).The minimum distance is thus found by a grid 

search with respect to all existing clusters, i = 1 1, and all feature vectors 

within the duration of an event I as defined in (20), 



Reference feature adaptation 

In order to track changes in the features of the different rhythms, 

adaptation of both £ imin (k) and i jm j^ (k) are desirable; here the event index k 
has been included for clarity. For £ jm j n (k), an exponentially updated average is 



dp (I) = (0(l) - Pi ) T S T 1 (<K0- Mi ) 



(24) 



25 



30 



used: 




(25) 
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where 

e(k)=<Knmin)-ML*.(k-1) ( 26 ) 



The exponential update factor y is confined to the interval 0 < y <1 . 
The inverse covariance matrix, £ jn ^ n (k), is estimated by exponential 

averaging of the new cluster difference matrix e(k) e T (k) using the update factor 
d-a), 



I. V) = 
• 



•mm 



v -1 



( a k U (k - 1)+ (1 - a>(k)s T (k) | (27) 
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Using the matrix inversion lemma 

A = B 1 + CD" 1 C T (28) 

A" 1 = B - BC(D + C T BC)- 1 C T B (29) 

15 and pairing the terms in (27) with the ones in (28), 

A = i imin (k) (30) 

B'aZin.inOc-l) (31) 

C - e(k) (32) 
20 D- 1 = (1 - a) (33) 

The inverse i jm jj] (k-1) in (27) may be computed without any matrix inversions as, 

(k-IMMs^kJa- 1 !- 1 (k-1) 

<i 44 row 'mm 

ST (k) = a~ 1 Sr 1 (k-1) : (34) 

'min 



(1-ar 1 +e T (k)a"Ti (k-1)e(k) 

'min 



25 By utilizing the matrix inversion lemma, the computational complexity 

of the operation is reduced from 0(P 3 ) to 0(P 2 ). 
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Algorithm initialization 

In leader-follower clustering, clusters are initialized as they become 
needed. This feature is convenient since it does not automatically introduce any 
unused clusters. However, it also puts demands on the algorithm to be able to 
5 create new clusters when necessary, and also to terminate clusters either not 
used for long or with only a few events. Initially, the total number of clusters, I, is 
equal to one. The algorithm is initialized by assigning the parameter vector 
0(n y ), which maximizes the test statistic in (14) of the first event, to the first 
cluster, cf. (23), 

10 

A,=<t>(n,) (35) 



15 



20 



25 



For the general case, <t>(n k ) is composed of a subset of the representative 
functions together with the preceding RR-interval, 



<D(n k ) = 



9s(n k ) 



L At k . 



(36) 



Where § s (ic k ) is a subset of the most discriminating elements in 0 s (A k ) . 

Note that an event is defined by its depolarization wave. Consequently, At 
is not included in the arrival time estimation o , but instead computed afterwards. 
The time continuous notation of At* is preferred since it results in a suitable 
magnitude similar to the normalized morphological information in 6 S . 

The inverse correlation matrix estimate is initialized in the same way for all 
clusters; a simple solution is to set it equal to a scaling of the identity matrix I, 



-1 _ 



= 81 



(37) 



where 5 is a design parameter. The complete clustering algorithm for 
organized events is presented in the flow chart in Fig. 2. 
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Reduced-complexity clustering 

In order to develop a more efficient algorithm in terms of performance 
versus power consumption and to evaluate the power consumption itself, a simple 




measure of the computational complexity can be used, namely, the total number of 
multiplications. This quantity represents a much more complex operation than do 
additions. In this algorithm, where most operations are of the nature "multiply- 
accumulate", the number of additions is of the same order as the number of 
5 multiplications and may thus be neglected without significant loss of accuracy. 

In order to reduce the computational complexity, focus is put on reducing 
the number of multiplications. The dominant contributions of multiplication 
operations are found in (24) and (34) which require P(P + 1)and P/2 (3P + 5) 

multiplications, respectively, considering certain symmetry properties of . 

10 Furthermore, one division is required in (34). However, according to (20), (24) is 
performed IK times per event while (34) is performed only once per event. 

Based on the above performance figures, a few approximations can 
be identified: 

o to use only the peak(s) in T(x(n)) instead of a complete grid search, 
15 o to use a simplified I"? , and 

o to use a likelihood based search sequence over the clusters, i.e., to 
start with the most likely cluster and to stop the search if a sufficiently 
small distance is found. 
Simplifying the grid search from spanning both samples and clusters 
20 in (20) to span over only clusters, 

d mm=mind?(n) (38) 
t 

the number of multiplications may be reduced by a factor K to IP(P 
+1 ). Due to sensitivity in T(x( n)) , the feature vectors resulting in the peaks for two 
different events may differ significantly, this simplification is likely to result in more 
25 clusters. A useful compromise may instead be to use the coefficients from, e.g., 
the 3 largest peaks in T(x(n)) resulting in 3IP(P + 1) multiplications per event. 

Since a cardiac event lasts for longer time than one sample those 
samples which give the filter coefficients which are most similar to the cluster 
reference are determined. A grid search over 40 msec is preferably made. In this 
30 way coefficients of greatest importance locally are determined. For this decision 
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T(x(n)) is considered and samples generating a peak are chosen, i.e. T(x(n-1)) < 
T(x(n)) > T(x(n+1 )), since they indicate the probability for the presence of an event. 

Another simplification, based on the approximation of orthogonality 

between the different wavelet scales as well as the RR interval, is to simplify £[~ 1 

S by only including its diagonal elements in the adaptation, 

£r 1 **" 1 £r 1 + ('-«)diag(s(k) E T (k)) (39) 

'min rrvn 

where diag(A) returns the diagonal elements of the square matrix A. By using this 
approximation, the number of multiplications used for the estimation of £j~ 1 is 

reduced to 3P. Additionally, the distance computation in (20) is simplified and may 
10 be reduced to 2IKP multiplications per event. 

A reasonable assumption is often that successive events originate 

from the same rhythms.Considering this knowledge, it would be sufficient to 

compute the distances for the previously selected cluster. In doing so, the number 

of multiplications in (38) are reduced even further. 
15 Table 1 presents the different detector versions as defined by their 

distinguishing featuresand shows computational complexity for the different 

versions of the clustering algorithm. 



Features 



Version 



-1 

imtn 



Search alg. 



Complexity C 



A 
B 

C 
D 



Full 
Diagonal 

Full 
Diagonal 



Interval 
Interval 

3 peak 
3 peak 



IKP(P + 1)+~ (3P + 5) 

2IKP+ 3P 
3IP(P + 1)+~ (3P + 5) 
6IP+ 3P 
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The total computational complexity, C, as reflected by the number of 
multiplications for the different algorithm versions, is presented in Fig. 3 as a 
function of I. In figure 3 the solid line shows an algorithm version A, dashed line a 




version B, dotted line a version C, and a version D by dash-dotted line), for a 
feature vector with (a), P = 4, and (b), P = 7 elements. 

Results 

5 The results are obtained by studying the performance of the following 

quantities: 

o algorithm versions, as presented in Table 1, and 
o noise tolerance, for noise-free signals and for signals with 
background noise of SO^V RMS. 
10 The following parameter settings have been used (unless otherwise 

stated), 

a = 1.05 y = 0.025 6 = 50 K = 40 (49) 

15 It is noted that a" 1 is chosen to offer faster adaption than y. The reason 

for that is that the initial estimate ir 1 is likely to be less accurate than the initial 
estimate pj . Also, 6 is chosen to have the same order of magnitude as the steady 

state eigenvalues . 

It should be noted that the different algorithm versions are not fully 
20 comparable in terms of performance for a specific due to the differences in 

distance computation. In versions B and D, where a diagonal t." 1 is used, the 

'min 

lack of non-diagonal information results in a nonorthogonal distance which is 
larger than the orthogonal one. Since versions C and D make use of a limited 
search, the minimum distance found may differ from the global minimum 
25 distance for the event. Both these algorithmic differences imply an increase in 
clustering quality for a certain value of n, however, at the expense of more 
introduced clusters. 

Evaluation measures 

30 The two main quantities evaluated are clustering performance and 

computational complexity, see Fig. 3; these two quantities are in general 
conflicting. In the evaluation, a cluster is assigned to that cardiac rhythm which 
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contains the most events in the cluster. This rhythm is denoted as the dominant 
rhythm within the cluster With respect to the dominant rhythm, the cluster is 
defined as a correct cluster, whereas it is erroneous to all other rhythms. If a 
rhythm is found to be dominant in more than one cluster, such clusters are first 

5 merged in the performance evaluation. The number of events in the correct cluster 
which belongs to the i:th dominant rhythm is equal to No(i), while the number of 
events belonging to any other false rhythm in the cluster is equal to N F (i). The 
number of events of the dominant rhythm which are not classified in a correct 
cluster, i.e., missing, is equal to N M (i). 

10 The performance of the algorithm is evaluated in terms of probability 

of a correct clustering of an event, PccO). and probability of a dominant event in a 
cluster, Pde(0. The first parameter is expressed as the share of correctly clustered 
events within a rhythm and is, using the above parameters, defined by 

Pcc(i)= Np(0 (41) 
W N D (i) + N M (i) 

15 The second parameter may be expressed as the share of dominant 

rhythm events within a cluster, and is defined as 

^E(0= M (42) 
N d (i) + N f (i) 

For the case when a rhythm completely lacks a correct cluster, Pde(0 
is undefined; the rhythm is then excluded from subsequent statistical 

20 computations. Averaging Pcc(i) andP D E(0 over all clusters results in the global 
performance measures Pec and Pde, respectively. 

It should be pointed out that Pcc(i) and PdeO) reach their maximal 
value of 1 when q is sufficiently small such that the number of clusters equals the 
total number of beats evaluated. This is, of course, a highly undesirable solution 

25 although performance, as expressed in (41) and (42), will be excellent. For this 
reason, the total number of clusters, I, is a crucial parameter to be considered. 
Here, the average T is used together with the minimum and maximum number of 
clusters for a case, l min , and l m ax, respectively. 

The power consumption of the algorithm is an important parameter 

30 which determines the pacemaker life span. In this study, power consumption is 
approximated by the computational complexity defined above as the total number 
of multiplications. As shown, the computational complexity depends mainly on 




three parameters: the feature vector length, P, the search interval, K and the 
cluster size, I. 

Noise-free signal clustering performance 

5 Figure 4 presents clustering performance in terms of Pde and Pec , 

see Fig. 4 (a), and 1mm, T and l ma * . see Fig. 4 (b). Thus in figure 4(a) clustering 
performance is shown as depending on T for Pec (dark bars) and Pde (bright bars) 
for noise-free signals. In figure 4(b), the spread between the different cases is 
shown as, from left to right, l m m, T and W The presented algorithm versions are 

10 found in Table 1 and three values of have been chosen for comparison; 3, 4 and 5. 
Versions A and B perform similarly for all three cases, and achieve Pde = 1 and 
p cc = 1 for T = 4 and T = 5, respectively, by creating an acceptable number of 
clusters. However, it can be seen from Fig. 4 (b) that a large difference in the 
number of initialized clusters between different cases is present. For version B, a 

15 slight increase in both T and l max is observed for T = 4. However, this increase 

is more due to an unfortunate step-like behaviour in the results for the given T than 
for any significant decrease in performance. The values of r\ used in Fig. 4 for the 
different versions are shown in Table 2. 

Table 2: Values of n used in Fig. 4. 

Algorithm version 



A 


B 


C 


D 


3.6 


4.2 


5.4 


5.8 


4.6 


4.8 


7.0 


7.2 


9.6 


9.6 


12.0 


10.8 
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Versions C and D perform slightly worse than do versions A and B. 
Contrary to the latter ones, neither version achieves Pde = 1 or P C c = 1 for the 
presented T, see figure 4(a). 

Figure 5 (a) presents the clustering performance for version A in detail 
25 and its dependence on q. The most noteworthy result is that both P DE = 1 and Pcc 
= 1 for q < 7. It is also clear that clustering performance deteriorates rapidly for n 
>10. The increase in P C c for very high values of q is due to that not all rhythms are 
allocated to a dominant cluster and are therefore disregarded in the performance 
computation. 
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In Fig. 5 (b), l min , T and Imaxare presented as a function of the 
clustering threshold. Thus figure 5(a) shows clustering performance in terms of 
Pde (solid line) and Pec (dashed line) as a function of n for noise-free signals, and 
figure 5(b) the corresponding U n , T and Imax- For n < 3, the number of initiated 
5 clusters increases rapidly with decreasing q. Within the interval 3 < n <7, all 
rhythms initiate at least one cluster, while for n >7, this is not true for all cases. 
Removing the \worst case", the above is true for the other cases for q <10. 

The clustering performance is exemplified for one set of parameters in 
Fig .6 using n = 7. Correct clustering with minimal number of clusters is achieved 

10 for two cases while an extra cluster is initialized for two cases. The reason for the 
extra cluster is that, for case 2, the temporal search interval is chosen too small for 
the third morphology from the left resulting in an extra cluster. For case 5 an actual 
difference in morphologies both on the up and down slopes of the dominant peak 
is discernible between the first and fourth clusters from the left. 

15 Clustering for a concatenated electrogram is presented for case 3 in 

Fig. 7, where three distinct rhythm classes can be discerned, viz. normal sinus 
rhythm followed by supraventricular tachycardia and atrial flutter. The EGM is 
shown together with the clusters, represented by o, x and +, respectively, assigned 
for each event.. The different rhythms result in three clusters. 

20 Figure 8 is a flow chart illustrating the function in broad outline of an 

embodiment of the apparatus according to the invention. At block 2 cardiac event 
features are extracted in the form of wavelet coefficients, and the event is 
detected, at 4, 6. At block 8 is checked whether the detected event is member of a 
labelled cluster. If so, the event is added to a class, at 10, and actions associated 

25 with that class are performed, at 12. 

If the event is not member of a labelled class it is checked if it is a 
member of an existing cluster, at 14. tf so, the event the event is added to a class, 
at 16, and it is checked if it is possible to label the cluster, 18, and if so the cluster 
is labeled, at 20. 

30 If the event is not a member of an existing cluster, block 14, a new 

cluster is created as described above fitted to the detected cardiac event, at 22. 

Thus according to the invention clustering events in the EGM is 
performed for use in implantable CRM devices, like heart stimulators. The 
invention is based on feature extraction in the wavelet domain whereupon the 




features are clustered based on the Mahalanobis distance criterion. According to 
advantageous embodiments of the apparatus according to the invention 
simplifications of the technique is proposed in order to reduce computational 
complexity to obtain a more implementationally feasible solution. 

5 If the apparatus according to the invention is to be used for longer 

periods of data analysis, large clusters, although old, may be desirable to be kept 
in some way, while the oldest cluster may be selected to be removed if the 
application is based on shorter time frames. Also, due to the short data lengths, 
testing of such algorithms would be of limited value. 

10 By combining detector/clusterer with labelling rules of a classifier a 

complete detector/classifier is obtained with the possibility to more accurate 
therapy. 

The labelling need not be done in real time and probably more than 
one event will be needed to label a cluster. Once the cluster has been labelled 

15 using clinical terms, the actions associated with the particular class will be carried 
out immediately, i.e. in real time. Thus, the rules needed to label the cluster are 
not used in identifying the event itself. 

The rules used to label the clusters are based on characteristics of the 
different possible events. Instead of the exact rules, the characteristics are 

20 consequently described. _ 

Figure 9 is a blockdiagram of an embodiment of an implantable heart 
stimulator provided with the apparatus according to the invention. Electrodes 30, 
32 implanted in the heart 34 of a patient are connected by a lead 36 to an A/D 
converter 40, via a switch 38 serving as overvoltage protection for the A/D 

25 converter 40. In the A/D converter 40 the signal is A/D converted and the digital 
signal is supplied to a wavelet detector 42. 

The detector 42 decides whether a cardiac event is present or not as 
described earlier. Wavelet coefficients are calculated as well. Parameters of the 
detector 42 are programmable from the stimulator microprocessor 44. At the 

30 detection the coefficients and the RR information are forwarded to the clusterer 46 
in which it is determined to which cluster the detected cardiac event belongs, as 
described previously. The clusterer 46 is preferably of a leader-follower type and 
also the cluster parameters are programmable from the microprocessor 44. 




By the microprocessor 44 suitable therapy is decided depending on 
assigned cluster for the detected event and possible a priori knowledge about 
arrhythmia associated with the cluster in question. Thus it is possible to distinguish 
e.g. ventricular tachycardia from a sinus tachycardia by comparing the parameters 

5 with a known normal sinus rhythm. Parameters of the sinus tachycardia are then 
supposed to be similar to those of a normal sinus rhythm, whereas parameters of 
ventricular tachycardia differ significantly. 

As an alternative the decision rules can be trained from a number of 
rhythms and the resulting rules are then used on test data, see Weichao Xu et al„ 

1 0 "New Bayesian Discriminator for Detection of Atrial Tachyarrhythmias", 

DOI:10.1161/01.CIR.0000012349.14270.54, pp.1472-1479, January 2002. It is 
then possible to decide that a certain position of the cluster indicates e.g. a sinus 
rhythm, etc. Also this technique can be based on analysis of the feature vector for 
a cluster, and it is possible to decide if the beat is broad or narrow, large or small, 

15 or if the rhythm is regular or irregular. 

The stimulator in figure 9 also includes a pulse unit 48 with associated 
battery 50 for delivery of stimulation pulses to the patient's heart 34 depending on 
the clustering evaluation. 

The implantable stimulator shown in figure 9 includes telemetry means 52 
20 with antenna 54 for communication with external equipment, like a programmer. 
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CLAIMS 

1 . An apparatus for analysing cardiac events detected in electrograms, 
EGMs, characterized in that a feature extraction means is provided to derive 

5 features of said cardiac events for discriminating different kinds of detected 
cardiac events, and in that a clusterer is provided to group cardiac events with 
similar features into a cluster, defined by predetermined cluster features. 

2. The apparatus according to claim 1 , characterized in that said feature 
10 extraction means is adapted to determine a feature vector describing waveform 

characteristics of cardiac event EGM signals. 

3. The apparatus according to claims 1 or 2, characterized in that said 
feature extraction means is adapted to perform morphological analysis on cardiac 

1 5 event EGM signals for the discrimination of different kinds of detected cardiac 
events. 

4. The apparatus according to any of the preceding claims, characterized In 
that said feature extraction means is adapted to determine the cardiac rate as one 

20 of said features. 

5. The apparatus according to claim 2, characterized in that said clusterer is 
adapted to determine the distance between said feature vector and corresponding 
cluster feature vectors to assign the cardiac event in question to that cluster which 

25 results in a minimum distance. 

6. The apparatus according to any of the preceding claims, each cluster 
being defined by a cluster center p; and a covariance matrix Ij for the respective 
cluster features, characterized in that said clusterer is adapted to determine a 

30 distance function d } 2 between event feature vector and said cluster center pi. 



7. The apparatus according to claim 6, characterized in that said clusterer is 
adapted to calculate said distance by using Mahalanobis distance criterion. 
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8. The apparatus according to claim 6 or 7, characterized in that said 
clusterer is adapted to determine the minimum of said distance by a grid search 
over the duration of the event. 

5 9. The apparatus according to any of the preceding claims, characterized in 
that an integrating means is provided to integrate said distance over a predeter- 
mined period of time. 

1 0. The apparatus according to any of the preceding claims, characterized in 
10 that said clusterer is adapted to update cluster features according to a predeter- 
mined rule on the basis of the difference between cluster features and event 
features. 

1 1 . The apparatus according to any of the preceding claims, characterized in 
15 that said clusterer is adapted to create a new cluster, if the difference between 

features of a detected cardiac event and features of existing clusters exceeds a 
predetermined limit, by setting features of a new cluster equal to the event 
features of the cardiac event in question. 

20 1 2. The apparatus according to any of the preceding claims, characterized in 
that said clusterer is adapted to terminate clusters not, or only rarely, being 
provided with detected cardiac events within a predetermined time period. 

13. The apparatus according to any of the preceding claims, characterized in 
25 that said clusterer is adapted to perform a likelihood based search sequence over 

the clusters to determine the minimum of said distance. 

14. The apparatus according to any of the claims 8-12, characterized in that 
said clusterer is adapted to to span said grid search only over the clusters. 

30 

15. The apparatus according to any of the claims 1 - 14, characterized in that 
said clusterer is adapted to calculate the distance of a considered cardiac event 
from the previously selected cluster. 
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16. The apparatus according to any of the preceding claims, characterized In 
that a classifier is provided to associate the clusters with different specific cardiac 
rhythms according to predetermined rules. 

17. A heart stimulator, characterized by an apparatus according to any of the 
preceding claims for on-line arrhythmia detection and control means for controlling 
the therapeutic stimulation depending on said arrhythmia detection. 
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ABSTRACT 

An apparatus for analysing cardiac events detected in electrograms, 
EGMs. comprises a feature extraction means (42) provided to derive features of 
said cardiac events for discriminating different kinds of detected cardiac events. A 
5 clusterer (46) is provided to group cardiac events with similar features into a 
cluster, defined by predetermined cluster features. A heart stimulator is provided 
with such an apparatus for arrhythmia detection and control means (44) for 
controlling the therapy depending on the arrhythmia detection. 



10 

(Figure 9) 
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